% Agent-based Model (ABM) code used to model the 2012 Republican's
% invisible primary by Nwokora & Brouzet (Nwokora, Zim, and Davy Brouzet. 
% "The Invisible Primary in an Agent-Based Model: Ideology, Strategy, and 
% Competitive Dynamics." Political Research Quarterly (2022):
% 10659129221107567. )
% Please read the license agreement before use


%% Formatting options
clearvars;close all;clc;
set(groot,'defaultFigureColor','w');
set(groot,'defaultLineLineWidth',2.5);
set(groot,'defaultAxesLineWidth',2);
set(groot,'defaultAxesFontSize',26);
set(groot,'defaultAxesTickLabelInterpreter','latex');
set(groot,'defaultAxesTickLength',[0.02 0.025]);
set(groot,'defaultAxesBox','on');
set(groot,'defaultFigureUnits','Normalized');
set(groot,'defaultFigurePosition',[0 0 0.5 0.7]);
set(groot,'DefaultFigureColormap',colormap(jet(64)));
set(groot,'DefaultColorbarTickLabelInterpreter','latex');
set(groot,'defaultTextInterpreter','latex');
set(groot,'defaultScatterLineWidth',2.5);
set(groot,'defaultScatterMarker','o');
set(groot,'defaultScatterMarkerEdgeColor','k');
set(groot,'defaultScatterMarkerFaceColor','k');
set(groot,'defaultScatterSizeData',100);
set(groot,'defaultContourLineWidth',2);
color_vec = get(groot,'DefaultAxesColorOrder');
color_vec(8,:) = [0,0,0];
set(groot,'DefaultAxesColorOrder',color_vec);
colormap(jet);
SizeFont = 30;
close all

%% Inputs to the ABM

n_ensemble = 2500;          %Number of times to the run ABM        

%Voters
SD_voters = 0.33;           %Standard deviations of voters
space_size_1 = 2;           %First dimension size
space_size_2 = 2;           %Second dimension size
correlation = 0.5;          %Correlation of voters distribution in economic and cultural dimensions
n_voters = 1000;            %Number of voters

%Candidates     
phi_vec = [1,0,0;0,1,0;0,0,1];                            %Strategy (1: Sticker, 2:Aggregator, 3: Hunter)
phi = [3,3,1,2,3,3,3,3];                                  %Strategy of each candidate
phi_1 = phi_vec(1,phi);                                   %Proportion of candidate being Sticker
phi_2 = phi_vec(2,phi);                                   %Proportion of candidate being Aggregator
phi_3 = 1.0 - phi_1 - phi_2;                              %Proportion of candidate being Hunter
gamma = SD_voters*[0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1];      %Displacement for every candidate (default is 10% of voters' SD)
kappa = [1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0];                %Satisficying level
delta_val = 11;                                           %Transition period to apply initial valence change (2 months)
val_start = [1.0,1.0,1.0,1.0,1.0,1.0,0.5,0.5];            %Valence starting value when candidate enters race (0.5)
val_1 = [1,1,1,1,1,1,1,1/1.4];                            %Baseline valence
delta_val_change = 4;                                     %Number of weeks for which the valence is changed
val_change = [1,8,8,7,2,2,0];                             %Candidate that has the valence advantage for each delta_val period
val_2 = [0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5];                %Valence change
n_it = length(val_change);


%Figures
show_graph='no';           %Show graphs
show_num='yes';            %Show number of supporters of graph
create_mov='no';           %Create movie

%% Voters' distribution, candidates' position and polls

%Get empirical voter's distributions
load('Density_Voters.mat','CF_vec','PL_vec','Density_CF','Density_PL');

%Order of candidates: Bachman, Gingrich, Huntsman, Paul, Romney, Santorum, Cain, Perry
name_abbr = ['B';'G';'H';'P';'R';'S';'C';'P'];
name = ['Bachman ';'Gingrich';'Huntsman';'Paul    ';'Romney  ';'Santorum';'Cain    ';'Perry   '];
n_parties = length(name);

%Give initial position of candidates
CF_init=[1.158;0.976;0.327;1.606;0.617;0.835;1.140;0.879];
PL_init=[0.895;0.908;0.065;0.714;0.261;0.999;0.827;0.676]*space_size_2;
pos_init=horzcat(CF_init,PL_init);

%Polls results
t_start=ones(n_parties,1);
t_start(7)=round(5);                 %Delayed start for Cain (5th week)
t_start(8)=round(15);                %Delayed start for Perry (15th week)
t1=datetime(2011,03,14,0,0,0);
t2=datetime(2012,01,16,0,0,0);
date_polls=t1:caldays(7):t2;
time_it=length(date_polls);
polls_tmp=csvread('Data_polls.csv',1,1,[1,1,length(date_polls),9]); %Read data                    
polls_norm=polls_tmp(:,1:8)*100./repmat(polls_tmp(:,9),1,8);
ENC_polls=1./sum((polls_norm/100).^2,2);                            %Compute ENC
init_results=polls_norm(1,1:n_parties)';

%% Perform simulation

%----------------------- Run n_ensemble indepedent ABM runs -------------------
path=zeros(2,n_parties,time_it,n_ensemble);
dist=zeros(n_parties,time_it,n_ensemble);
psize=zeros(n_parties,time_it,n_ensemble);
init_error=zeros(n_parties,n_ensemble);
Rp=zeros(time_it,n_ensemble);
disp('Running ABM...');
for n=1:n_ensemble
    
    clc;
    disp([num2str(n/n_ensemble*100,3) '%']);
    %Call Laver's algorithm
    [parties,voters]=Politics_Laver_2D(n_parties,name,name_abbr,phi_1,phi_2,phi_3,...
                                       gamma,kappa,delta_val,delta_val_change,val_start,val_1,val_change,val_2,n_it,t_start,pos_init, ...
                                       n_voters,Density_CF,Density_PL,CF_vec,PL_vec,space_size_1,space_size_2,correlation, ...
                                       time_it,polls_norm,show_num,show_graph,create_mov,color_vec,n);

    %Post-processing        
    for i=1:n_parties
       path(1,i,:,n)=parties.xcor(i,1:time_it);
       path(2,i,:,n)=parties.ycor(i,1:time_it);
       dist(i,:,n)=sqrt( parties.xcor(i,1:time_it).^2 +...
       parties.ycor(i,1:time_it).^2 )/SD_voters;
       psize(i,:,n)=parties.psize(i,:);
       init_error(i,n)=abs(parties.psize(i,1)/n_voters*100-init_results(i));
    end
    Rp(:,n) = voters.Rp;

end

%% Further post-processing and showing results

%Compute mean absolute difference MAD
MAD=zeros(n_parties,1);
for i=1:n_parties
   MAD(i)=mean(abs(polls_norm(t_start(i):time_it,i)'-mean(psize(i,t_start(i):time_it,:),3)/n_voters*100));
end
fitness = sqrt(mean(MAD.^2));
fitness_candidate = MAD;
disp(['Averaged MAD (mean absolute difference): ' num2str(sqrt(mean(MAD.^2)))]);
for i=1:n_parties
   disp([name(i,:) ' ' num2str(MAD(i))]);
end

%Show time traces
figure(11)
for i=1:n_parties
 plot(date_polls(1:time_it),mean(psize(i,:,:),3)/n_voters*100);hold on;
 if i==9
   plot(date_polls(1:time_it),squeeze(psize(i,:,:)/n_voters*100),'Color',color_vec(i,:),'Linewidth',0.25);hold on;
 end
end
set(gca,'YLim',[0 50]);
xlabel('Time','Interpreter','latex');
ylabel('Parties size (\%)','Interpreter','latex');
legend(name);

%Show time traces with error bars and polls
figure(12)
for i=1:n_parties
  subplot(3,3,i)
  plot(date_polls(1:time_it),polls_norm(1:time_it,i),'Color','k','LineStyle','-');hold on;
  plot(date_polls(1:time_it),mean(psize(i,:,:),3)/n_voters*100,'LineStyle','-','Color','r');hold on;
  for k=1:time_it-1
     %Shows the standard deviation, can have some issues plotting depeding
     %on your MATLAB version and/or installed options
     %patch([date_polls(k) date_polls(k+1) date_polls(k+1) date_polls(k)],[mean(psize(i,k,:),3)-std(psize(i,k,:),[],3) mean(psize(i,k+1,:),3)-std(psize(i,k+1,:),[],3) ...
     %       mean(psize(i,k+1,:),3)+std(psize(i,k+1,:),[],3) mean(psize(i,k,:),3)+std(psize(i,k,:),[],3)]/n_voters*100,'r','EdgeColor','r','FaceAlpha',0.25,'LineStyle','none')  
  end
  set(gca,'YLim',[0 50]);
  title(name(i,:));
end
set(gcf,'Position',[0.0 0.0 1.0 1.0]);

%Representativeness
figure(13)
plot(date_polls(1:time_it),mean(Rp,2));hold on;
xlabel('Time','Interpreter','latex');
ylabel('Representativeness','Interpreter','latex');
set(gcf,'Position',[0.0 0.0 0.35 0.45]);
ylim([-0.1 0]);

%Effective number of candidates
figure(14)
ENC_sim=1./sum((mean(psize(:,:,:),3)/n_voters).^2,1);
plot(date_polls(1:time_it),ENC_polls,'Color','k','LineStyle','-');hold on;
plot(date_polls(1:time_it),ENC_sim,'LineStyle','--');hold on;
xlabel('Time','Interpreter','latex');
ylabel('ENC','Interpreter','latex');
set(gcf,'Position',[0.0 0.525 0.35 0.45]);

return



